%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Wave packet propagation integrating the Schrödinger equation,
%  1D time-splitting spectral method.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;								Constant;
figure(1);				        	clf;

a  =-10e-6;							b  =10e-6;
M  =512;							dz =(b-a)/M;								%grid points
j  =[0:M-1];						z  =a+j*dz;
l  =[-M/2:M/2-1];					k0 =2*pi*l/(b-a);
dt =5e-6;							wz =2*pi*50;								%time step
									V  =87*uA/2*wz^2*z.^2;						%trap
%epl=100e6;							V  =hP*2*pi*1.5e3*(1+atan((z-10e-6)*epl)/pi-atan((z+10e-6)*epl)/pi);
%lmd=5e-6;							V  =hP*2*pi*1.5e4*exp((z-20e-6)/lmd);
E  =hP*wz*(0+.5);
k  =.0*sqrt(2*87*uA*E)/hP;			aHO=sqrt(hP/87/uA/wz);
Psi=exp(-(z/aHO).^2/2).*exp(i*k*z);	Psi=Psi/sqrt(sum(abs(Psi).^2)*dz);	

skp=5;
for kk=1:1000
%	[kk,sqrt(sum(abs(Psi).^2)*dz)]
	Psi_x =exp(-i/hP*V*dt/2) .* Psi;
	Psi_xf= fft(Psi_x);
	Psi_y =exp(i*hP/2/87/uA*k0.^2*dt) .* Psi_xf;
	Psi_yf=ifft(Psi_y);
	Psi   =exp(-i/hP*V*dt/2) .* Psi_yf;
	if((kk/skp)==floor(kk/skp))
		subplot(121),plot(z*1e6,V/hP/wz,'b',z*1e6,E/hP/wz,'m',z*1e6,abs(Psi).^2*1e-5+E/hP/wz,'r')
		set(gca,'xLim',[min(z*1e6) max(z*1e6)],'yLim',[0 5],'FontSize',16,'LineWidth',1.4)%
		xlabel('\fontsize{18}{\itx}  (\mum)')
		ylabel('\fontsize{18}{\itV/h}  (kHz)')
		pause(.01);
	end
end